Estimating the impact of non-pharmaceutical interventions on COVID-19 infection spread through mobility data

We aim at answering the question :

How does is the implementation of different strategies affecting the rates of COVID-19 infection ?

Summary

Through the Covid-19 pandemic, many countries have adopted mobility restriction policies, such as lockdown. Major technology companies have released mobility datasets, that describe the mobility change, relative to a baseline, across different categories of places or modes of transportation.

We observe that countries cluster based on the mobility data, revealing patterns and similarities between government measures across different countries. We build a model that predicts mobility change based on non-pharmaceutical interventions (NPIs).
Additionally, we build a bayesian model, that leverages an epidemiological compartmental model, in order to estimate the real effective reproduction number $R_t$ through time as a function of mobility data.

This opens the possibility to build an end-to-end pipeline from NPIs to $R_t$ and compartments, e.g. number of hospitalized, critical and deceased individuals, in order to generate what-if scenarios, useful to evaluate governement policies both restrospectively and prospectively.

process

Datasets

Measures dataset

We leverage the Coronavirus Government Response Tracker dataset from Oxford Blavatnik School of Government.

We are interested in the Closure and containment indicators, which are described in the following section. Note that when no measure is in place, the indicator is at zero.

Measures description

C1 School closing
Record closings of schools and universities

  • 1: recommend closing
  • 2: Require closing (only some levels or categories, eg just high school, or just public schools)
  • 3: Require closing all levelsNo data

C2 Workplace closing
Record closings of workplaces

  • 1: recommend closing (or work from home)
  • 2: require closing (or work from home)for some sectors or categories of workers
  • 3: require closing (or work from home) all-but-essential workplaces (eg grocery stores, doctors)

C3 Cancel public events
Record cancelling public events

  • 1: Recommend cancelling
  • 2: Require cancelling

C4 Restrictions on gatherings
Record the cut-off size for bans on private gatherings

  • 1: Restrictions on very large gatherings (the limit is above 1000 people)
  • 2: Restrictions on gatherings between 100-1000 people
  • 3: Restrictions on gatherings between 10-100 people
  • 4: Restrictions on gatherings of less than 10 people

C5 Close public transport
Record closing of public transport

  • 1: Recommend closing (or significantly reduce volume/route/means of transport available)
  • 2: Require closing (or prohibit most citizens from using it)

C6 Stay at home requirements
Record orders to “shelter-in-place” and otherwise confine to home

  • 1: recommend not leaving house
  • 2: require not leaving house with exceptions for daily exercise, grocery shopping, and ‘essential’ trips
  • 3: Require not leaving house with minimal exceptions (e.g. allowed to leave only once every few days, or only one personcan leave at a time, etc.)

C7 Restrictions on internal movement
Record restrictions on internal movement

  • 1: Recommend closing(or significantly reduce volume/route/means of transport)
  • 2: Require closing (or prohibit most people from using it)

C8 International travel controls
Record restrictions on international travel

  • 1: Screening
  • 2: Quarantine arrivals from high-risk regions
  • 3: Ban on high-risk region

todo: insérer ici des viz du dataset oxford

Mobility datasets

We gathered data from Google and Apple, describing respectively 6 and 3 mobility categories.

Google

  • Retail and recreation
    • restaurants, cafes, shopping centers, theme parks, museums, libraries, and movie theaters.
  • Grocery and pharmacy
    • grocery markets, food warehouses, farmers markets, specialty food shops, drug stores, and pharmacies.
  • Parks
    • national parks, public beaches, marinas, dog parks, plazas,and public gardens
  • Transit stations
    • public transport hubs such as subway, bus, and train stations
  • Workplaces
  • Residential

The baseline is the median value, for the corresponding day of the week, during the 5-week period Jan 3–Feb 6, 2020.

Apple

  • Driving
  • Walking
  • Transit

relative volume of directions requests per country/region, sub-region or city compared to a baseline volume on 13 January 2020

todo: insérer ici des viz de mobilité

Relationship between measures and mobility

Data fetching
/home/horace/miniconda3/envs/uncoverun/lib/python3.8/site-packages/IPython/core/interactiveshell.py:3254: DtypeWarning: Columns (3) have mixed types.Specify dtype option on import or set low_memory=False.
  if (await self.run_code(code, result,  async_=asy)):
retail_and_recreation grocery_and_pharmacy parks transit_stations workplaces residential driving transit walking
iso_code date
ABW 2020-02-15 0.02 0.03 0.12 0.26 -0.02 -0.02 NaN NaN NaN
2020-02-16 -0.01 -0.09 0.04 0.01 -0.01 0.00 NaN NaN NaN
2020-02-17 -0.01 -0.01 0.16 -0.04 -0.20 0.06 NaN NaN NaN
2020-02-18 -0.02 -0.05 0.08 0.02 -0.13 0.04 NaN NaN NaN
2020-02-19 0.01 -0.08 0.03 0.04 -0.14 0.04 NaN NaN NaN
... ... ... ... ... ... ... ... ... ... ...
ZAF 2020-05-04 NaN NaN NaN NaN NaN NaN -0.5624 NaN -0.5293
2020-05-05 NaN NaN NaN NaN NaN NaN -0.5412 NaN -0.5008
2020-05-06 NaN NaN NaN NaN NaN NaN -0.5258 NaN -0.4535
2020-05-07 NaN NaN NaN NaN NaN NaN -0.5423 NaN -0.4638
2020-05-08 NaN NaN NaN NaN NaN NaN -0.5483 NaN -0.4635

13299 rows × 9 columns

Mobility, a good proxy for measures

The intuition we had, working on this problem, was that mobility could be a useful proxy for measures. In order to verify this hypothesis we decided to take a few approaches.

  • First we compared the mobility timeseries to see if countries with similar policies have similar mobilities
  • Then we learn the relationship between measures and mobility through bayesian inference.
Countries used for learning,  ['DNK', 'NLD', 'JPN', 'CZE', 'BEL', 'PHL', 'SVK', 'TWN', 'CAN', 'ESP', 'USA', 'CHE', 'GBR', 'BRA', 'NOR', 'NZL', 'FIN', 'DEU', 'AUS', 'MEX', 'EST', 'FRA', 'IRL', 'LUX', 'SWE']
<matplotlib.axes._subplots.AxesSubplot at 0x7f62d69cc130>

Mobility analysis

  • Let's cluster the different countries by typology of mobility
0it [00:00, ?it/s]
dict_keys(['GUM', 'JOR', 'PRT', 'PER', 'YEM', 'DNK', 'ECU', 'AFG', 'HRV', 'UKR', 'ITA', 'IDN', 'GRC', 'LAO', 'NLD', 'MNG', 'CIV', 'MMR', 'ALB', 'THA', 'MUS', 'KEN', 'CMR', 'LBY', 'JPN', 'TZA', 'MLI', 'CZE', 'BEL', 'SVN', 'PHL', 'SVK', 'KOR', 'TWN', 'CAN', 'PRI', 'HUN', 'SGP', 'ESP', 'BEN', 'URY', 'BGD', 'NIC', 'USA', 'CHE', 'NGA', 'ABW', 'ARE', 'POL', 'PNG', 'KAZ', 'CRI', 'AUT', 'KWT', 'ZAF', 'BRB', 'HKG', 'SLV', 'CPV', 'PAK', 'GBR', 'BRA', 'PAN', 'TTO', 'KGZ', 'ROU', 'NOR', 'GHA', 'MOZ', 'NZL', 'CHL', 'RWA', 'MAR', 'COL', 'IRQ', 'ISL', 'BFA', 'MYS', 'PRY', 'VEN', 'BGR', 'ZMB', 'LKA', 'NER', 'NAM', 'BLZ', 'FIN', 'QAT', 'RUS', 'DEU', 'DOM', 'MDA', 'VNM', 'AUS', 'ISR', 'AGO', 'BOL', 'BIH', 'ZWE', 'BHR', 'IND', 'GAB', 'SAU', 'OMN', 'ARG', 'MEX', 'HND', 'LBN', 'EST', 'SRB', 'FRA', 'GTM', 'TUR', 'IRL', 'JAM', 'LUX', 'BWA', 'EGY', 'SWE', 'UGA'])
27it [03:18,  7.37s/it]
<matplotlib.axes._subplots.AxesSubplot at 0x7f62d6e90f10>

Similarity Analysis

We see some interesting strong similarities between countries. Like Italy and France that have been hit hard by Covid-19 and enforced a total lockdown. But how can we visualize these typologies better?

/home/horace/miniconda3/envs/uncoverun/lib/python3.8/site-packages/seaborn/matrix.py:619: ClusterWarning: scipy.cluster: The symmetric non-negative hollow observation matrix looks suspiciously like an uncondensed distance matrix
  linkage = hierarchy.linkage(self.array, method=self.method,
<seaborn.matrix.ClusterGrid at 0x7f62d73cb550>

Cluster Analysis

With this new point of view, we are able to detect countries that have indeed taken similar measures leading to similar mobility. For instance this approach shows that Japan and Sweden have similar mobility's profiles. They have both avoided strong lockdown. On the other hand Italy, France and Great Britain belong to the same cluster of drastic mobility loss. Finally we can also detect countries that have decided to handle the crisis with intermediate measures, like Denmark, Germany or Norway.Therefore it seems that this clusterization validates the use of mobility as a convenient proxy for measures.

Bayesian regression

Let $M_{t, i} \geq 0$ be here the relative change in mobility, expressed here differently to ease modelling.
In this section, we thus have $M_{t_i} = 1$ when there is no change in mobility, and $M_{t, i} = 0$ when a 100% diminution has occured (theoretic).

For a mobility category $i$, time $t$, we model $M_{t, i} = \mathbb{E}[m_{t, i}]$ with:

$$M_{t, i} \sim \mathcal{N}(m_{t, i}, \sigma)$$

where the prior for $\sigma$ is $\sigma \sim \mathrm{Exponential}(1) $.

We model:

$$m_{t, i} = \exp(- \sum_k \alpha_{k, i} I_{k, t})$$

with

  • $I_{k, t} \in \mathbb{N}$ ordinal indicator of the measure $k$ at time $t$
  • $\alpha_{k, i} \geq 0$ to be inferred, with prior $\alpha_{k, i} \sim \mathrm{Gamma}(\frac{1}{N}, 1) - \frac{\log 1.05}{N}$ where $N$ is the number of measures

For simplicity, we do not include the residential category, which has opposite correlation with mobility. We also remove the parks series, that is quite noisy in some countries.

We fit this model on multiple countries with the same parameters, in order to extract a country-independent effect of measures on mobility.

array(['ABW', 'AFG', 'AGO', 'ALB', 'AND', 'ARE', 'ARG', 'AUS', 'AUT',
       'AZE', 'BDI', 'BEL', 'BEN', 'BFA', 'BGD', 'BGR', 'BHR', 'BIH',
       'BLZ', 'BMU', 'BOL', 'BRA', 'BRB', 'BRN', 'BWA', 'CAN', 'CHE',
       'CHL', 'CHN', 'CIV', 'CMR', 'COD', 'COG', 'COL', 'CPV', 'CRI',
       'CUB', 'CYP', 'CZE', 'DEU', 'DJI', 'DMA', 'DNK', 'DOM', 'DZA',
       'ECU', 'EGY', 'ESP', 'EST', 'ETH', 'FIN', 'FRA', 'GAB', 'GBR',
       'GHA', 'GMB', 'GRC', 'GRL', 'GTM', 'GUM', 'GUY', 'HKG', 'HND',
       'HRV', 'HUN', 'IDN', 'IND', 'IRL', 'IRN', 'IRQ', 'ISL', 'ISR',
       'ITA', 'JAM', 'JOR', 'JPN', 'KAZ', 'KEN', 'KGZ', 'KOR', 'KWT',
       'LAO', 'LBN', 'LBR', 'LBY', 'LKA', 'LSO', 'LUX', 'MAC', 'MAR',
       'MDA', 'MDG', 'MEX', 'MLI', 'MMR', 'MNG', 'MOZ', 'MRT', 'MUS',
       'MWI', 'MYS', 'NAM', 'NER', 'NGA', 'NIC', 'NLD', 'NOR', 'NZL',
       'OMN', 'PAK', 'PAN', 'PER', 'PHL', 'PNG', 'POL', 'PRI', 'PRT',
       'PRY', 'PSE', 'QAT', 'RKS', 'ROU', 'RUS', 'RWA', 'SAU', 'SDN',
       'SGP', 'SLE', 'SLV', 'SMR', 'SOM', 'SRB', 'SSD', 'SUR', 'SVK',
       'SVN', 'SWE', 'SWZ', 'SYC', 'SYR', 'TCD', 'THA', 'TKM', 'TTO',
       'TUN', 'TUR', 'TWN', 'TZA', 'UGA', 'UKR', 'URY', 'USA', 'UZB',
       'VEN', 'VNM', 'YEM', 'ZAF', 'ZMB', 'ZWE'], dtype=object)
sample: 100%|██████████| 1000/1000 [02:24<00:00,  6.90it/s, 1023 steps of size 4.76e-03. acc. prob=0.93]
                mean       std    median      5.0%     95.0%     n_eff     r_hat
         a      0.01      0.00      0.01      0.01      0.02    265.13      1.02
alpha[0,0]      0.27      0.03      0.28      0.22      0.32     23.05      1.06
alpha[0,1]      0.15      0.03      0.15      0.11      0.19    144.88      1.00
alpha[0,2]      0.06      0.01      0.06      0.04      0.09     86.67      1.01
alpha[0,3]      0.00      0.01      0.00      0.00      0.01     57.78      1.02
alpha[0,4]      0.07      0.01      0.07      0.05      0.09    549.61      1.00
alpha[0,5]      0.11      0.01      0.11      0.09      0.12    759.41      1.00
alpha[0,6]      0.14      0.01      0.14      0.12      0.17    200.95      1.00
alpha[1,0]      0.29      0.03      0.29      0.23      0.34     30.08      1.04
alpha[1,1]      0.08      0.04      0.08      0.00      0.12     87.20      1.01
alpha[1,2]      0.00      0.00      0.00      0.00      0.00    428.85      1.00
alpha[1,3]      0.00      0.00      0.00      0.00      0.00    411.54      1.00
alpha[1,4]      0.10      0.02      0.10      0.06      0.13    127.77      1.01
alpha[1,5]      0.14      0.01      0.14      0.12      0.16    756.12      1.00
alpha[1,6]      0.18      0.02      0.18      0.15      0.20    328.06      1.00
alpha[2,0]      0.10      0.02      0.09      0.06      0.12    146.24      1.02
alpha[2,1]      0.00      0.01      0.00      0.00      0.00    406.86      1.02
alpha[2,2]      0.01      0.02      0.00      0.00      0.04     11.10      1.07
alpha[2,3]      0.00      0.00      0.00      0.00      0.00    503.98      1.00
alpha[2,4]      0.00      0.00      0.00      0.00      0.00    345.85      1.00
alpha[2,5]      0.00      0.00      0.00      0.00      0.00    322.38      1.01
alpha[2,6]      0.00      0.00      0.00      0.00      0.00    190.55      1.00
alpha[3,0]      0.01      0.01      0.00      0.00      0.03      7.52      1.17
alpha[3,1]      0.00      0.00      0.00      0.00      0.00    269.63      1.02
alpha[3,2]      0.00      0.00      0.00      0.00      0.00    504.04      1.00
alpha[3,3]      0.00      0.00      0.00      0.00      0.00    329.66      1.00
alpha[3,4]      0.00      0.00      0.00      0.00      0.00    484.02      1.00
alpha[3,5]      0.00      0.00      0.00      0.00      0.00    222.73      1.00
alpha[3,6]      0.02      0.01      0.02      0.01      0.03     67.95      1.03
alpha[4,0]      0.00      0.00      0.00      0.00      0.00    398.67      1.00
alpha[4,1]      0.00      0.00      0.00      0.00      0.00    421.61      1.00
alpha[4,2]      0.00      0.00      0.00      0.00      0.00    398.78      1.00
alpha[4,3]      0.04      0.01      0.04      0.02      0.06    301.25      1.00
alpha[4,4]      0.00      0.00      0.00      0.00      0.00    473.50      1.00
alpha[4,5]      0.00      0.00      0.00      0.00      0.00    455.69      1.01
alpha[4,6]      0.00      0.00      0.00      0.00      0.00    417.40      1.00
alpha[5,0]      0.08      0.06      0.09      0.00      0.16     15.95      1.09
alpha[5,1]      0.36      0.05      0.36      0.29      0.45    213.74      1.00
alpha[5,2]      0.52      0.02      0.52      0.49      0.57    299.15      1.00
alpha[5,3]      0.26      0.01      0.26      0.25      0.28    206.68      1.00
alpha[5,4]      0.50      0.03      0.50      0.46      0.56    409.64      1.00
alpha[5,5]      0.13      0.02      0.13      0.10      0.16    764.02      1.00
alpha[5,6]      0.19      0.02      0.19      0.16      0.23    364.91      1.00
alpha[6,0]      0.00      0.00      0.00      0.00      0.00    450.85      1.00
alpha[6,1]      0.00      0.00      0.00      0.00      0.00    485.88      1.00
alpha[6,2]      0.00      0.00      0.00      0.00      0.00    501.90      1.00
alpha[6,3]      0.00      0.00      0.00      0.00      0.00    397.05      1.00
alpha[6,4]      0.00      0.01      0.00      0.00      0.00     55.21      1.01
alpha[6,5]      0.00      0.01      0.00      0.00      0.00    296.33      1.01
alpha[6,6]      0.00      0.00      0.00      0.00      0.00    495.24      1.00
alpha[7,0]      0.00      0.00      0.00      0.00      0.00    506.78      1.00
alpha[7,1]      0.00      0.00      0.00      0.00      0.00    391.32      1.00
alpha[7,2]      0.00      0.00      0.00      0.00      0.00    481.28      1.00
alpha[7,3]      0.02      0.00      0.02      0.01      0.02    199.10      1.00
alpha[7,4]      0.02      0.00      0.02      0.01      0.03    497.14      1.00
alpha[7,5]      0.02      0.00      0.02      0.02      0.03    511.51      1.01
alpha[7,6]      0.03      0.00      0.03      0.03      0.04    535.30      1.00
  sigma[0]      0.12      0.00      0.12      0.12      0.13    464.72      1.01
  sigma[1]      0.19      0.01      0.19      0.18      0.20    597.04      1.00
  sigma[2]      0.14      0.00      0.14      0.14      0.15    690.22      1.00
  sigma[3]      0.08      0.00      0.08      0.07      0.08    617.86      1.01
  sigma[4]      0.11      0.00      0.11      0.10      0.11    999.10      1.00
  sigma[5]      0.07      0.00      0.07      0.07      0.08    519.16      1.00
  sigma[6]      0.08      0.00      0.08      0.07      0.08    624.02      1.00

Number of divergences: 2
<ipython-input-160-dbd6fcd3c4ea>:15: RuntimeWarning: More than 20 figures have been opened. Figures created through the pyplot interface (`matplotlib.pyplot.figure`) are retained until explicitly closed and may consume too much memory. (To control this warning, see the rcParam `figure.max_open_warning`).
  fig, ax = plt.subplots(figsize=(8, 5))
<ipython-input-160-dbd6fcd3c4ea>:15: RuntimeWarning: More than 20 figures have been opened. Figures created through the pyplot interface (`matplotlib.pyplot.figure`) are retained until explicitly closed and may consume too much memory. (To control this warning, see the rcParam `figure.max_open_warning`).
  fig, ax = plt.subplots(figsize=(8, 5))
<ipython-input-160-dbd6fcd3c4ea>:15: RuntimeWarning: More than 20 figures have been opened. Figures created through the pyplot interface (`matplotlib.pyplot.figure`) are retained until explicitly closed and may consume too much memory. (To control this warning, see the rcParam `figure.max_open_warning`).
  fig, ax = plt.subplots(figsize=(8, 5))
<ipython-input-160-dbd6fcd3c4ea>:15: RuntimeWarning: More than 20 figures have been opened. Figures created through the pyplot interface (`matplotlib.pyplot.figure`) are retained until explicitly closed and may consume too much memory. (To control this warning, see the rcParam `figure.max_open_warning`).
  fig, ax = plt.subplots(figsize=(8, 5))
<ipython-input-160-dbd6fcd3c4ea>:15: RuntimeWarning: More than 20 figures have been opened. Figures created through the pyplot interface (`matplotlib.pyplot.figure`) are retained until explicitly closed and may consume too much memory. (To control this warning, see the rcParam `figure.max_open_warning`).
  fig, ax = plt.subplots(figsize=(8, 5))
<ipython-input-160-dbd6fcd3c4ea>:15: RuntimeWarning: More than 20 figures have been opened. Figures created through the pyplot interface (`matplotlib.pyplot.figure`) are retained until explicitly closed and may consume too much memory. (To control this warning, see the rcParam `figure.max_open_warning`).
  fig, ax = plt.subplots(figsize=(8, 5))
<ipython-input-160-dbd6fcd3c4ea>:15: RuntimeWarning: More than 20 figures have been opened. Figures created through the pyplot interface (`matplotlib.pyplot.figure`) are retained until explicitly closed and may consume too much memory. (To control this warning, see the rcParam `figure.max_open_warning`).
  fig, ax = plt.subplots(figsize=(8, 5))
<ipython-input-160-dbd6fcd3c4ea>:15: RuntimeWarning: More than 20 figures have been opened. Figures created through the pyplot interface (`matplotlib.pyplot.figure`) are retained until explicitly closed and may consume too much memory. (To control this warning, see the rcParam `figure.max_open_warning`).
  fig, ax = plt.subplots(figsize=(8, 5))
<ipython-input-160-dbd6fcd3c4ea>:15: RuntimeWarning: More than 20 figures have been opened. Figures created through the pyplot interface (`matplotlib.pyplot.figure`) are retained until explicitly closed and may consume too much memory. (To control this warning, see the rcParam `figure.max_open_warning`).
  fig, ax = plt.subplots(figsize=(8, 5))
<ipython-input-160-dbd6fcd3c4ea>:15: RuntimeWarning: More than 20 figures have been opened. Figures created through the pyplot interface (`matplotlib.pyplot.figure`) are retained until explicitly closed and may consume too much memory. (To control this warning, see the rcParam `figure.max_open_warning`).
  fig, ax = plt.subplots(figsize=(8, 5))
<ipython-input-160-dbd6fcd3c4ea>:15: RuntimeWarning: More than 20 figures have been opened. Figures created through the pyplot interface (`matplotlib.pyplot.figure`) are retained until explicitly closed and may consume too much memory. (To control this warning, see the rcParam `figure.max_open_warning`).
  fig, ax = plt.subplots(figsize=(8, 5))
<ipython-input-160-dbd6fcd3c4ea>:15: RuntimeWarning: More than 20 figures have been opened. Figures created through the pyplot interface (`matplotlib.pyplot.figure`) are retained until explicitly closed and may consume too much memory. (To control this warning, see the rcParam `figure.max_open_warning`).
  fig, ax = plt.subplots(figsize=(8, 5))
<ipython-input-160-dbd6fcd3c4ea>:15: RuntimeWarning: More than 20 figures have been opened. Figures created through the pyplot interface (`matplotlib.pyplot.figure`) are retained until explicitly closed and may consume too much memory. (To control this warning, see the rcParam `figure.max_open_warning`).
  fig, ax = plt.subplots(figsize=(8, 5))
<ipython-input-160-dbd6fcd3c4ea>:15: RuntimeWarning: More than 20 figures have been opened. Figures created through the pyplot interface (`matplotlib.pyplot.figure`) are retained until explicitly closed and may consume too much memory. (To control this warning, see the rcParam `figure.max_open_warning`).
  fig, ax = plt.subplots(figsize=(8, 5))
<ipython-input-160-dbd6fcd3c4ea>:15: RuntimeWarning: More than 20 figures have been opened. Figures created through the pyplot interface (`matplotlib.pyplot.figure`) are retained until explicitly closed and may consume too much memory. (To control this warning, see the rcParam `figure.max_open_warning`).
  fig, ax = plt.subplots(figsize=(8, 5))
<ipython-input-160-dbd6fcd3c4ea>:15: RuntimeWarning: More than 20 figures have been opened. Figures created through the pyplot interface (`matplotlib.pyplot.figure`) are retained until explicitly closed and may consume too much memory. (To control this warning, see the rcParam `figure.max_open_warning`).
  fig, ax = plt.subplots(figsize=(8, 5))
<ipython-input-160-dbd6fcd3c4ea>:15: RuntimeWarning: More than 20 figures have been opened. Figures created through the pyplot interface (`matplotlib.pyplot.figure`) are retained until explicitly closed and may consume too much memory. (To control this warning, see the rcParam `figure.max_open_warning`).
  fig, ax = plt.subplots(figsize=(8, 5))
<ipython-input-160-dbd6fcd3c4ea>:15: RuntimeWarning: More than 20 figures have been opened. Figures created through the pyplot interface (`matplotlib.pyplot.figure`) are retained until explicitly closed and may consume too much memory. (To control this warning, see the rcParam `figure.max_open_warning`).
  fig, ax = plt.subplots(figsize=(8, 5))
<ipython-input-160-dbd6fcd3c4ea>:15: RuntimeWarning: More than 20 figures have been opened. Figures created through the pyplot interface (`matplotlib.pyplot.figure`) are retained until explicitly closed and may consume too much memory. (To control this warning, see the rcParam `figure.max_open_warning`).
  fig, ax = plt.subplots(figsize=(8, 5))
<ipython-input-160-dbd6fcd3c4ea>:15: RuntimeWarning: More than 20 figures have been opened. Figures created through the pyplot interface (`matplotlib.pyplot.figure`) are retained until explicitly closed and may consume too much memory. (To control this warning, see the rcParam `figure.max_open_warning`).
  fig, ax = plt.subplots(figsize=(8, 5))
<ipython-input-160-dbd6fcd3c4ea>:15: RuntimeWarning: More than 20 figures have been opened. Figures created through the pyplot interface (`matplotlib.pyplot.figure`) are retained until explicitly closed and may consume too much memory. (To control this warning, see the rcParam `figure.max_open_warning`).
  fig, ax = plt.subplots(figsize=(8, 5))
<ipython-input-160-dbd6fcd3c4ea>:15: RuntimeWarning: More than 20 figures have been opened. Figures created through the pyplot interface (`matplotlib.pyplot.figure`) are retained until explicitly closed and may consume too much memory. (To control this warning, see the rcParam `figure.max_open_warning`).
  fig, ax = plt.subplots(figsize=(8, 5))

Models

SEIR-HCD Model

Model built on compartmental models. Each of the following letters represents a compartment of the population of a country:

S - susceptible
E - exposed
I - infectious
R - recovered
H - hospitalized
C - critical
D - deceased

We normalize these count numbers by the total population of the country.

\begin{align} \frac{dS}{dt} &= - \frac{R_t}{T_{inf}} I S \\ \frac{dE}{dt} &= \frac{R_t}{T_{inf}} I S - \frac{E}{T_{inc}} \\ \frac{dI}{dt} &= \frac{E}{T_{inc}} - \frac{I}{T_{inf}} \\ \frac{dR}{dt} &= m_a \frac{I}{T_{inf}} + (1 - c_a)\frac{H}{T_{hosp}} \\ \frac{dH}{dt} &= (1 - m_a) \frac{I}{T_{inf}} + (1 - f_a)\frac{C}{T_{crit}} - \frac{H}{T_{hosp}} \\ \frac{dC}{dt} &= c_a \frac{H}{T_{hosp}} - \frac{C}{T_{crit}} \\ \frac{dD}{dt} &= f_a \frac{C}{T_{crit}} \end{align}

Parameters used in the model

$R_t$ = reproduction number at time t.
Typical 3.6* at t=0

Transition times

T_inc = average incubation period. Typical 5.6* days
T_inf = average infectious period. Typical 2.9 days
T_hosp = average time a patient is in hospital before either recovering or becoming critical. Typical 4 days
T_crit = average time a patient is in a critical state (either recover or die). Typical 14 days

Fractions
These constants are likely to be age specific (hence the subscript a):

m_a = fraction of infections that are asymptomatic or mild. Assumed 80% (i.e. 20% severe)
c_a = fraction of severe cases that turn critical. Assumed 10%
f_a = fraction of critical cases that are fatal. Assumed 30%

*Averages taken from https://www.kaggle.com/covid-19-contributions

Priors for these parameters

  • The Conversation article:

    • Incubation: 5 days
    • Infectious period: 2 days in incubation, 6 days from symptoms onwards
    • Illness: 10 days
    • NB: For COVID-19, there is emerging evidence to suggest the infectious period may start 1 to 3 days before you develop symptoms
      • We can reduce incubation time by 2 days and increase infectious time by 2
    • article30196-1/fulltext)
      • The peak viral load of patients with MERS-CoV and SARS-CoV infections occurs at around 7–10 days after symptom onset
      • The median interval between symptom onset and hospitalisation was 4 days (range 0–13)
  • Imperial

    • T_inc: Gamma(5.1, 0.86)
      • The infection-to-onset distribution is Gamma distributed with mean 5.1 days and coefficient of variation 0.86.
    • T_inf + T_hosp + T_crit: Gamma(17.8, 0.45)
      • The onset-to-death distribution is also Gamma distributed with a mean of 17.8 days and a coefficient of variation 0.45.
    • T_inc + T_inf + T_hosp + T_crit
      • The infection-to-death distribution is therefore given by: π∼Gamma(5.1,0.86) + Gamma(17.8,0.45)

Reproduction number

let $\mathrm{m}_{t, i}$ be the reduction of mobility in the category $i$, relatively to a baseline (before the pandemic).
Hence $\mathrm{m}_{t, i} > -1$.

We adopt the strong hypothesis that each mobility category presents the same transmission dynamics. Let $R_{t, i}$ be defined by the following linear relationship: $$ R_{t, i} = (1 + \mathrm{m}_{t, i})R_0 - \mathrm{m}_{t, i}R_1 $$ for any time $t$ through the pandemic, with $R_0$, $R_1$ to be estimated.

We model the effective reproduction number with a weighted mean for each mobility category:

$$R_t = \sum_i \alpha_i R_{t, i}$$

with $\alpha_i$ to be estimated, where

  • $\sum_i \alpha_i = 1$
  • $\forall i, \alpha_i > 0$

We choose the following prior : $\forall i, \alpha_i \sim \mathrm{Gamma}(1, .5)$

Effective reproduction number and death modelling

Parameterization

We use the following parameterizations:

  • (concentration, rate) for the $\mathrm{GammaPoisson}$ distribution.
  • (mean, std) for the $\mathrm{Gamma}$ distribution
  • (mean, sample size) for the $\mathrm{Beta}$ distribution

Target

Let $D_t$ be the number of death at time $t$ for a given country. We model $d_t = \mathbb{E}[D_t]$

We sample $$D_t \sim \mathrm{GammaPoisson}(\psi, \frac{\psi}{d_t})$$ where $\psi \sim \mathcal{N}^+(0, 5)$

Compartment initialization

We seed the model with the following, assuming $N$ is the country population:

  • $I_0 \sim \mathrm{Gamma}(100, 50) / N$
  • $S_0 = 1 - I_0$
  • $E_0, R_0, H_0, C_0, D_0 = 0$

Other priors

$R_0 \sim \mathcal{N}^+(3.28, \kappa_0)$

$R_1 \sim \mathcal{N}^+(0.7, \kappa_1)$

where $\kappa_0, \kappa_1 \sim \mathcal{N}^+(0, .5)$

$T_{inc} \sim \mathrm{Gamma}(5.6, .86)$

$T_{inf} \sim \mathrm{Gamma}(2.9, 1)$

$T_{hosp} \sim \mathrm{Gamma}(4, 1)$

$T_{crit} \sim \mathrm{Gamma}(14, 1)$

$m_a \sim \mathrm{Beta}(0.8, \phi_m)$

$c_a \sim \mathrm{Beta}(0.1, \phi_c)$

$f_a \sim \mathrm{Beta}(0.35, \phi_f)$

where $\phi_m, \phi_c, \phi_f \sim \mathrm{Gamma}(7, 2)$

Inference settings

We use a multi-country setting where we train the model on a pool of countries and then predict the number of deceased individuals, as well as $R_t$, on another country.

In order to do this, the model has to learn country-independent parameters, since the only varying factor in countries is the mobility.

Data fetching

/home/horace/miniconda3/envs/uncoverun/lib/python3.8/site-packages/IPython/core/interactiveshell.py:3254: DtypeWarning: Columns (3) have mixed types.Specify dtype option on import or set low_memory=False.
  if (await self.run_code(code, result,  async_=asy)):
train countries: ['DEU', 'ESP', 'ITA']
test countries: ['DNK', 'FRA', 'SWE']

Inference

sample: 100%|██████████| 400/400 [56:05<00:00,  8.41s/it, 1023 steps of size 5.69e-06. acc. prob=0.83] 
                     mean       std    median      5.0%     95.0%     n_eff     r_hat
       alpha[0]      2.80      0.00      2.80      2.80      2.81      3.04      1.81
       alpha[1]      2.25      0.00      2.25      2.25      2.25      4.56      1.20
       alpha[2]      0.23      0.00      0.23      0.23      0.23      3.34      2.35
       alpha[3]      1.37      0.00      1.37      1.37      1.37      2.79      2.61
       alpha[4]      2.46      0.00      2.46      2.46      2.47      8.44      1.28
       alpha[5]      5.74      0.00      5.74      5.74      5.74      4.47      1.24
       alpha[6]      0.30      0.00      0.30      0.30      0.30      5.11      1.04
            c_a      0.24      0.00      0.24      0.24      0.24      4.21      1.57
            f_a      0.68      0.00      0.68      0.68      0.68      2.64      2.35
         i_init      3.46      0.00      3.46      3.46      3.46      9.67      1.00
         kappa0      1.35      0.00      1.35      1.35      1.35      6.00      1.01
         kappa1      2.81      0.00      2.81      2.81      2.81      5.46      1.01
            m_a      0.32      0.00      0.32      0.32      0.32      5.13      1.36
            psi      7.76      0.01      7.76      7.74      7.78      2.63      2.32
             r0      1.87      0.00      1.87      1.87      1.87      5.83      1.33
             r1      4.94      0.00      4.94      4.93      4.94      3.52      1.32
  sample_size_c      0.52      0.00      0.52      0.52      0.52      4.14      1.25
  sample_size_f      0.19      0.00      0.19      0.19      0.19      6.79      1.01
  sample_size_m      2.48      0.00      2.48      2.48      2.48      3.06      2.22
         t_crit      0.38      0.00      0.38      0.38      0.38      4.84      1.43
         t_hosp      0.15      0.00      0.15      0.15      0.15      2.87      2.06
          t_inc      0.83      0.00      0.83      0.83      0.83      2.90      2.01
          t_inf      3.31      0.00      3.31      3.31      3.31      3.60      1.69

Number of divergences: 0

Results

Populations of all compartments over time

<ipython-input-55-5b700d381edc>:36: RuntimeWarning: More than 20 figures have been opened. Figures created through the pyplot interface (`matplotlib.pyplot.figure`) are retained until explicitly closed and may consume too much memory. (To control this warning, see the rcParam `figure.max_open_warning`).
  plt.subplots()
<ipython-input-55-5b700d381edc>:43: RuntimeWarning: More than 20 figures have been opened. Figures created through the pyplot interface (`matplotlib.pyplot.figure`) are retained until explicitly closed and may consume too much memory. (To control this warning, see the rcParam `figure.max_open_warning`).
  plt.subplots()
<ipython-input-55-5b700d381edc>:43: RuntimeWarning: More than 20 figures have been opened. Figures created through the pyplot interface (`matplotlib.pyplot.figure`) are retained until explicitly closed and may consume too much memory. (To control this warning, see the rcParam `figure.max_open_warning`).
  plt.subplots()
<ipython-input-55-5b700d381edc>:43: RuntimeWarning: More than 20 figures have been opened. Figures created through the pyplot interface (`matplotlib.pyplot.figure`) are retained until explicitly closed and may consume too much memory. (To control this warning, see the rcParam `figure.max_open_warning`).
  plt.subplots()
<ipython-input-55-5b700d381edc>:43: RuntimeWarning: More than 20 figures have been opened. Figures created through the pyplot interface (`matplotlib.pyplot.figure`) are retained until explicitly closed and may consume too much memory. (To control this warning, see the rcParam `figure.max_open_warning`).
  plt.subplots()
<ipython-input-55-5b700d381edc>:36: RuntimeWarning: More than 20 figures have been opened. Figures created through the pyplot interface (`matplotlib.pyplot.figure`) are retained until explicitly closed and may consume too much memory. (To control this warning, see the rcParam `figure.max_open_warning`).
  plt.subplots()
<ipython-input-55-5b700d381edc>:36: RuntimeWarning: More than 20 figures have been opened. Figures created through the pyplot interface (`matplotlib.pyplot.figure`) are retained until explicitly closed and may consume too much memory. (To control this warning, see the rcParam `figure.max_open_warning`).
  plt.subplots()
<ipython-input-55-5b700d381edc>:36: RuntimeWarning: More than 20 figures have been opened. Figures created through the pyplot interface (`matplotlib.pyplot.figure`) are retained until explicitly closed and may consume too much memory. (To control this warning, see the rcParam `figure.max_open_warning`).
  plt.subplots()
<ipython-input-55-5b700d381edc>:43: RuntimeWarning: More than 20 figures have been opened. Figures created through the pyplot interface (`matplotlib.pyplot.figure`) are retained until explicitly closed and may consume too much memory. (To control this warning, see the rcParam `figure.max_open_warning`).
  plt.subplots()
<ipython-input-55-5b700d381edc>:43: RuntimeWarning: More than 20 figures have been opened. Figures created through the pyplot interface (`matplotlib.pyplot.figure`) are retained until explicitly closed and may consume too much memory. (To control this warning, see the rcParam `figure.max_open_warning`).
  plt.subplots()
<ipython-input-55-5b700d381edc>:43: RuntimeWarning: More than 20 figures have been opened. Figures created through the pyplot interface (`matplotlib.pyplot.figure`) are retained until explicitly closed and may consume too much memory. (To control this warning, see the rcParam `figure.max_open_warning`).
  plt.subplots()
<ipython-input-55-5b700d381edc>:43: RuntimeWarning: More than 20 figures have been opened. Figures created through the pyplot interface (`matplotlib.pyplot.figure`) are retained until explicitly closed and may consume too much memory. (To control this warning, see the rcParam `figure.max_open_warning`).
  plt.subplots()
<ipython-input-55-5b700d381edc>:36: RuntimeWarning: More than 20 figures have been opened. Figures created through the pyplot interface (`matplotlib.pyplot.figure`) are retained until explicitly closed and may consume too much memory. (To control this warning, see the rcParam `figure.max_open_warning`).
  plt.subplots()
<ipython-input-55-5b700d381edc>:36: RuntimeWarning: More than 20 figures have been opened. Figures created through the pyplot interface (`matplotlib.pyplot.figure`) are retained until explicitly closed and may consume too much memory. (To control this warning, see the rcParam `figure.max_open_warning`).
  plt.subplots()
<ipython-input-55-5b700d381edc>:36: RuntimeWarning: More than 20 figures have been opened. Figures created through the pyplot interface (`matplotlib.pyplot.figure`) are retained until explicitly closed and may consume too much memory. (To control this warning, see the rcParam `figure.max_open_warning`).
  plt.subplots()

With synthetic mobility data

todo

Conclusion

todo
NB: effectiveness of masks

Future work

  • We could leverage hospitalization and ICU data and feed it to the model, so that it learns based on observations of H, C and D compartments at once. This would reinforce the confidence one can have on the deceased compartment, and $R_t$, predictions, as well as enable the model to estimate other compartments more accurately. An example of dataset that could be used is the one from IHME that provides estimated numbers.

  • We could include other data sources as additional mobility measurements.

  • Combining to make scenarios: one could build a framework that has the following properties:
    • Input : measure, date of implementation, country
    • Output : mobility, $R_t$, deaths [+ infected, hospitalized, critical] time series

Here are some example questions of interest:

  • What if the lockdown had been decided one week later ? Two weeks earlier ?
  • What will happen if the lockdown is lifted in one week ? Will we see a spike in hospitals ?

These questions are of utmost importance for policy makers.

Sources

Data sources

  • Epidemiology
    • Our world in data: confirmed cases and deaths from ECDC (European center of Disease Control) + tests (collected by owid)
  • Mobility
    • Apple - 3 categories: walking, driving, transit
    • Google - 6 categories